%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
%tilde='/home/mwoolrich';
tilde='/Users/woolrich';
osldir=[tilde '/homedir/matlab/osl1.2'];    

addpath(osldir);
osl_startup(osldir);

spm eeg;
ca;

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

testdir='/Users/woolrich/homedir/matlab/osl_testdata_dir';

datadir=[testdir '/faces_subject1_data']; % directory where the data is

workingdir=[datadir]; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

% Set up the list of subjects and their structural scans for the analysis 
% Currently there is only 1 subject.
clear spm_files;

% set up a list of SPM MEEG object file names (we only have one here)
spm_files{1}=[workingdir '/spmfiles/espm8_meg1.mat'];

% set up a list of mris
structural_files{1}=[datadir '/structurals/struct1.nii'];

%%%%%%%%%%%%%%%%%%%
%% First do a TF analysis

S2=[];
S2.D_continuous=spm_files;
oat.source_recon.conditions={'Motorbike','Neutral face','Happy face','Fearful face'};
S2.freq_range=[1 40]; % frequency range in Hz
S2.time_range=[-0.2 0.3];
%S2.first_level_time_range=[0.1 0.22];

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
S2.design_matrix_summary=Xsummary;     

% contrasts to be calculated:
S2.contrast={};
S2.contrast{1}=[3 0 0 0]'; % motorbikes
S2.contrast{2}=[0 1 1 1]'; % faces
S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
S2.recon_method='beamform';
S2.gridstep=14; % in mm, using a lower resolution here than you would normally, for computational speed
S2.tf_method='hilbert';
S2.bc=[1 1 1];

oat = osl_setup_oat(S2);

%%%%%%%%%%%%%%%%%%%
%% RUN OAT
oat.source_recon.dirname=[workingdir '/subj1_results.oat'];
oat.first_level.cope_type='cope';
oat.first_level.bc=[1 1 1];

oat.to_do=[1 1 0];
oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%%%%%
%% OUTPUT SUBJECT'S NIFTII FILES
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[3]; % list of first level contrasts to output
S2.stats_dir='/Users/woolrich/homedir/matlab/faces_subject1_data/spmfiles/subj1_results.oat/subject1_first_level_dir';
[statsdir,times]=osl_save_nii_stats(S2);

%%%%%%%
%% do beamformer using ROI using envelope
mni_coords=[34,-70,-2];
gridstep=2;
oat.source_recon.mask_fname=osl_mnicoords2mnimask(mni_coords,gridstep,'hi');
oat.source_recon.time_range=[-0.2 0.8];
oat.first_level.freq_range=[2 20]; % frequency range in Hz
oat.first_level.time_range=[-0.15 0.6];
oat.first_level.baseline_timespan=[-0.15 -0.05];
oat.source_recon.dirname=[workingdir '/subj1_results_roi.oat'];
oat.first_level.tf_method='hilbert';
oat.first_level.num_freqs=30;
oat.first_level.bc=[1 1 0];
oat.first_level.cope_type='cope';

oat.first_level.space_average=1;
oat.to_do=[1 1 0];
oat = osl_run_oat(oat);

% load GLM result
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

% Plot the parameter estimates for our closest beamformed voxel over time, for 
% our two contrasts of interest:
figure;
con=3; imagesc(stats.times, stats.frequencies, squeeze((stats.cope(1,:,con,:)))');axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['cope' num2str(con)]);
figure;imagesc(stats.times, stats.frequencies, squeeze((stats.cope(1,:,con,:)./stats.stdcope(1,:,con,:)))');axis xy;
ylabel('frequency (Hz)'); xlabel('time (s)'); colorbar; title(['tstat' num2str(con)]);

%%%%%%%%%%%%%%%%%%%
%% SETUP BEAMFORMER AND FIRST-LEVEL GLM OAT connectivity using envelope
S2=[];
S2.D=spm_files;
S2.trigger={{'Motorbike'},{'Neutral face'},{'Happy face'},{'Fearful face'}};
S2.freq_range=[2 10]; % frequency range in Hz
S2.time_range=[-0.2 0.8];
%S2.first_level_time_range=[0.1 0.22];

% Xsummary is a parsimonious description of the design matrix.
% It contains values Xsummary{reg,cond}, where reg is a regressor no. and cond
% is a condition no. This will be used (by expanding the conditions over
% trials) to create the (num_regressors x num_trials) design matrix:
Xsummary={};
Xsummary{1}=[1 0 0 0];Xsummary{2}=[0 1 0 0];Xsummary{3}=[0 0 1 0];Xsummary{4}=[0 0 0 1];
S2.design_matrix_summary=Xsummary;

% contrasts to be calculated:
S2.contrast={};
S2.contrast{1}=[3 0 0 0]'; % motorbikes
S2.contrast{2}=[0 1 1 1]'; % faces
S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
S2.recon_method='beamform';
S2.gridstep=14; % in mm, using a lower resolution here than you would normally, for computational speed
%S2.connectivity_seed_mni_coord=[30,-56,-14];
S2.connectivity_seed_mni_coord=[34,-70,-2];
S2.connectivity_method=@osl_conn_hilbert_envelope_ds;
S2.first_level_freq_range=[2 10]; % frequency range in Hz
S2.first_level_time_range=[0.05 0.4];
S2.first_level_baseline_timespan=[-0.15 -0.05];
S2.oat_name=[workingdir '/subj1_results_conn_hilb_env.oat'];

oat = osl_setup_oat(S2);

%%%%%%%%%%%%%%%%%%%
%% RUN OAT connectivity using envelope

oat.first_level.cope_type='cope';
oat.first_level.connectivity_method=@osl_conn_hilbert_envelope;
oat.first_level.bc=[0 0 0]; % do not want to do baseline correction 
oat.first_level.time_moving_av_win_size=0.05;

oat.to_do=[0 1 0];
oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%%%%%
%% OUTPUT SUBJECT'S NIFTII FILES for connectivity using envelope
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1:3]; % list of first level contrasts to output
[statsdir,times]=osl_save_nii_stats(S2);

%%%%%%%%%%%%%%%%%
%% VIEW NIFTII RESULTS IN FSLVIEW
% We can now view the nifti images containing our GLM results in FSL, here
% we are runnin fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.

mni_2mm_brain=[getenv('FSLDIR') '/data/standard/MNI152_T1_2mm_brain']; 
        
% INSPECT THE RESULTS OF A CONTRAST IN FSLVIEW:
%S2.contrast{1}=[3 0 0 0]'; % motorbikes
%S2.contrast{2}=[0 1 1 1]'; % faces
%S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
con=3;
runcmd([ 'fslview ' mni_2mm_brain ' ' [statsdir '/cope' num2str(con) '_2mm'] ' ' [statsdir '/tstat' num2str(con) '_2mm']  ' &']);



%%%%%%%%%%%%%%%%%%%
%% RUN OAT using phase
oat.source_recon.dirname=[workingdir '/subj1_results.oat'];
oat.first_level.cope_type='cope';
oat.first_level.connectivity_method=@osl_conn_hilbert_plv;
oat.source_recon.freq_range=[11 17]; % frequency range in Hz
oat.source_recon.time_range=[-0.2 0.8];
oat.first_level.freq_range=[13 15]; % frequency range in Hz
oat.first_level.time_range=[0.05 0.15];
oat.first_level.baseline_timespan=[-0.1 0];

oat.first_level.bc=[1 1 1]; 

oat.to_do=[1 1 0];
oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%%%%%
%% OUTPUT SUBJECT'S NIFTII FILES
% Having run the GLM on our source space data, we would like to inspect the
% results for our single subject. 
% We can do this by saving the contrast of parameter estimates (COPEs) and 
% t-statistics for each of our contrasts to NIFTI images.

S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1,2,3]; % list of first level contrasts to output
S2.stats_dir='/Users/woolrich/homedir/matlab/faces_subject1_data/spmfiles/subj1_results.oat/subject1_first_level_dir';
[statsdir,times]=osl_save_nii_stats(S2);

%%%%%%%%%%%%%%%%%
%% VIEW NIFTII RESULTS IN FSLVIEW
% We can now view the nifti images containing our GLM results in FSL, here
% we are runnin fslview from the matlab command line, but you do not need
% to - you can run it from the UNIX command line instead.

mni_2mm_brain=[getenv('FSLDIR') '/data/standard/MNI152_T1_2mm_brain']; 
        
% INSPECT THE RESULTS OF A CONTRAST IN FSLVIEW:
%S2.contrast{1}=[3 0 0 0]'; % motorbikes
%S2.contrast{2}=[0 1 1 1]'; % faces
%S2.contrast{3}=[-3 1 1 1]'; % faces-motorbikes
con=1;
runcmd([ 'fslview ' mni_2mm_brain ' ' [statsdir '/cope' num2str(con) '_2mm'] ' ' [statsdir '/tstat' num2str(con) '_2mm']  ' &']);





runcmd([ 'fslview ' mni_2mm_brain ' ' [statsdir '/tstat' num2str(1) '_2mm'] ' ' [statsdir '/tstat' num2str(2) '_2mm'] ' ' [statsdir '/tstat' num2str(3) '_2mm'] ' &']);













  